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ABSTRACT 

In  this  report,  we  present  a  survey  of  problems  and 
solutions  in  the  area  of  target  tracking.  The  discussion  includes 
design  trade-offs,  performance  evaluation,  and  current  issues. 
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1. 


INTRODUCTION 


For  nearly  two  decades  the  target  tracking-tra jectory  esti¬ 
mation  problem  has  been  a  fruitful  applications  area  for  state 
estimation.  Many  problems  have  been  solved  yet  new  and  diversi¬ 
fied  applications  still  challenge  systems  engineers.  In  this 
report,  we  provide  a  review  of  the  subject  and  discuss  some  current 
issues. 

An  earlier  excellent  reference  on  tracking- trajectory  esti¬ 
mation  [1]  deals  with  tracking  algorithms  for  ballistic  re-entry 
vehicles.  Most  of  these  algorithms  are  also  applicable  to  a  more 
general  class  of  targets  such  as  tactical  missiles  and  airplanes. 
The  emphasis  of  this  paper  will  therefore  not  be  on  algorithmic 
details  but  rather  on  discussing  problem  areas  and  approaches. 
References  containing  estimation  algorithms  are  cited  for  more 
interested  readers.  Furthermore,  discussions  contained  in  this 
report  will  be  oriented  from  the  practitioner's  point  of  view  and 

consequently  rigorous  mathematical  terms  are  only  of  secondary 
interest. 

This  report  is  organized  as  follows.  The  fundamental  problem 
of  target  tracking,  approaches  and  some  design  tradeoffs  are 

reviewed  in  Section  2.  Four  approaches  for  tracking  targets, 
with  sudden  maneuvers  are  presented  in  Section  3.  The  discussion 
includes  a  comparison  of  their  relative  merits.  The  problem  of 
tracking  with  passive  sensors  (measurements  containing  only  line- 


of-sight  angles)  is  discussed  in  Section  4.  In  some  applications, 
one  is  confronted  with  large  scale  system  issues,  i.e.,  the 
existence  of  many  sensors  operating  in  a  multiple  target  environ¬ 
ment.  In  Section  5,  several  algorithms  for  processing  multiple 
sensor  data  and  an  algorithm  for  correlating  measurements  from 
multiple  sensors  are  presented.  The  problem  of  tracking  in  a 
multiple  target  environment  is  discussed  in  Section  6.  In  each 
tracking  application,  one  may  be  interested  in  performance  evalu¬ 
ation  without  resorting  to  Monte  Carlo  simulations.  Covariance 
analysis  techniques  are  outlined  in  Section  VII  for  this  purpose, 
namely,  the  polynomial  analysis,  the  Riccati  Equation,  and  the 
Cramer- Rao  Bound. 

Since  this  report  focuses  discussions  on  design  considerations, 
performance  trade-offs  and  approaches  to  given  problems,  several 
related  subjects  discussing  algorithmic  details  such  as  the  square 
root  [53]  or  U-D  factorization  algorithm  [54]  are  not  included. 

One  may  notice  that  certain  subjects  are  given  more  attention 
than  others,  this  is  due  primarily  to  personal  preferences.  In 
each  case  however,  we  try  to  provide  an  adequate  list  of  refer¬ 
ences  to  allow  a  indepth  study  for  interested  readers. 
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FUNDAMENTALS 


2.1.  Problem  Definition 

The  tracking  problem  is  a  state  estimation  problem,  i.e., 
assuming  the  state  of  a  target  evolves  in  time  according  to  the 
equation 

x  =  f  (x)  +  w  (2.1) 

and  the  corresponding  (discrete)  measurement  vector*  is  given  by 


Ik  =  -(-k*  +  -k  (2.2) 

where  w  and  v^  are  input  and  measurement  noise  processes,  respec¬ 
tively,  one  is  interested  in  estimating  the  target  states  based 
upon  all  measurements  z^,  £=l,...,k. 

We  make  these  following  remarks. 


(1)  Eq.  (2.1)  is  a  mathematical  model  representative  of 
the  target  dynamics.  The  state  vector  x,  usually 
contains  target  position,  velocity,  and  sometimes 
acceleration  as  state  variables.  In  some  sit¬ 
uations,  key  parameters  characterizing  important 
target  properties  are  also  included  as  state 
variables.  The  filter  designer  usually  has  the 
option  of  choosing  among  several  models  with 
different  level  of  complexity.  The  trade-off  is 
performance  versus  real-time  computational  require¬ 
ment. 

(2)  Eq.  (2.2)  is  the  measurement  equation  relating 
state  variables  to  measurement  variables.  When 

a  radar  is  used,  z.  has  at  least  three  components, 
i.e.,  range  and  two  angles.  If  a  passive  sensor 
(such  as  a  telescope)  is  used,  z .  only  contains 
two  angle  measurements.  K 

*Ihe  measurement  device  can  be  a  radar,  sonar,  telescope  (passive) ,  and  others. 

In  most  cases,  these  measurements  are  taken  in  discrete  times. 
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(3)  The  system  (input)  and  measurement  noise  processes 
w  and  Vfc,  respectively,  are  assumed  to  be  zero 
mean  white  noise  processes.  The  covariance  of 
w,  Q  is  selected  to  compensate  for  modeling 
errors  (discrepencies  between  (2.1)  and  the 
actual  process).  The  statistics  of  the  measure¬ 
ment  noise  process  vj^  should  also  be  selected  to 
represent  all  possible  excursions  such  as 
measurement  biases,  false  measurements,  etc. 


2.2  Basic  Approaches 


The  basic  tracking  filter  is  a  recursive  algorithm. 
During  time  tfc  to  tk+1,  the  state  estimate  is  computed  by 
integrating 

x  =  f(x)  (2.3) 


from  tk  to  tk+1 
new  measurement 
by 


using  xk^,k  as  the  initial  state.  At  time  tfc+1, 
£k+1  is  obtained,  the  state  estimate  is  updated 


a 


-k+l/k+1  "  — k+l/k  +  Kk+1 (— k+1  "  -(-k+l/k;' 


A  # 

where  x^^  is  the  estimate  of  x^  based  upon  measurements 
z k=l,...,j.  Two  questions  arise: 

(1)  How  does  one  choose  f(  )? 

(2)  How  is  the  filter  gain  Kk+1  computed? 

These  two  questions  which  appear  different  initially  are  actually 
intimately  related  and  are  discussed  below. 
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2.2.1  Target  Dynamics 

For  some  targets,  a  constant  velocity  (CV)  model  is 
sufficient,  i.e.,  the  state  vector  contains  six  variables 

x  =  [x1,x2,x3,x4,x5,x6]T  =  [x,x,y,y,z,z]T  (2.5) 

where  (x,y,x)  are  coordinates  of  a  Cartesian  coordinate  system 
used  to  describe  the  target  dynamics.  The  state  equations  are 


x.  =  x. , . 
l  l+l 


X. , .  =  w. 

1+1  1 


i=l,3,5 


(2.6) 


where  w^  is  a  process  noise  term  used  to  characterize  modeling 
errors. 

If  the  target  being  tracked  is  maneuvering  (accelerating) , 
a  constant  acceleration  (CA)  model  is  usually  used,  i.e.. 


•  ••  •  •  •  ••• 

x  =  [x,x,x,y,y,y ,z,z,z] 

The  state  equations  can  be  written  accordingly,  i.e.. 


(2.7) 


Xi  “  xi+l 


*i+i  -  *i+2' 


(2.8) 


x .  .  _  =  w. 

1+2  l 
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These  two  state  equations  are  also  referred  to  as  the  first  and 
second  order  polynomial  dynamics,  respectively. 


Note  also  that  equations  (2.6)  and  (2.8)  assume  complete 
decoupling  between  x,y,  and  z.  A  commonly  used  Cartesian  system 
has  the  coordinate  centered  at  the  sensor  location  with  the  x-axis 
pointing  east,  the  y-axis  pointing  north,  and  the  z-axis  per¬ 
pendicular  to  the  local  horizontal  plane. 

The  target  dynamic  is  sometimes  described  in  the  sensor 
coordinates.  For  example,  if  a  dish  radar  is  used  for  tracking, 

the  measurement  variables  include  range  (r) ,  azimuth  (a) ,  and 
elevation  (e) .  A  constant  velocity  target  dynamic  model  decoupled 
in  r,  a,  and  e  results  in  the  following  equations. 


r  =  w 


•  • 

a  =  w_ 


(2.9) 


Where  wr,  wq,  and  wg  are  process  noise  terms  representing  modeling 
errors  in  r,  a,  and  e  directions,  respectively.  The  advantage 
of  using  (2.9)  is  that  one  is  required  to  construct  three  two- 
dimensional  filters  instead  of  one  six-dimensional  filter.  Con¬ 
siderable  computational  savings  result.  This  model  however,  is 
inconsistent  with  the  assumption  the  target  dynamics  is  decoupled 


in  the  Cartesian  coordinate,  eqs.  (2.6),  and  (2.8).  A  more 
appropriate  model  is  to  retain  the  coupling  among  r,  a  and  e. 

This  results  in  the  following  set  of  equations. 

•  •  .*2  *2  2 

r  =  r(e  +  a  cos  e) 

a  =  -2— a  +  2  a  e  tan  e  (2.10) 

r 

.  .2 

•  •  r*  3. 

e  =  -2— e  -  —  sin  2e 
r  1 

We  note  that  using  measurement  variables  as  state 
variables  may  result  in  more  accurate  state  estimates  because 
this  makes  the  measurement  equations  linear  (see  for  example,  the 
discussion  of  [1])«  The  Cartesian  coordinates  however,  appeal 
intuitively  and  provide  easier  interpretation  of  target  motion. 

These  models  are  often  used  for  tracking  airplanes  and 
tactical  missiles  and  are  also  used  for  tracking  ballistic  missiles. 
In  some  re-entry  (RV)  vehcile  applications,  however ,  the  vehicle  aero¬ 
dynamic  parameters  such  as  the  ballistic  and  lifting  coefficients 
must  be  estimated  in  real  time.  This  requires  a  set  of  nonlinear 
differential  equations  to  describe  RV  motion.  A  simplified*  model 
is 

a(t)  =  ^pv2(t)(jUd  +  atufc  +  ac^cl  +  9  (2.11) 

where  a(t)  is  the  total  acceleration  applied  on  the  vehicle 

♦This  is  a  simplified  model  neglecting  Coriolis  and  Centrifugal 
forces.  They  are  usually  included  using  the  vector  sum  as  in  (2.11). 
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p 


is  the  air  density 


v(t)  is  the  magnitude  of  the  vehicle  velocity 

3  is  the  ballistic  coefficient 

u  is  the  unit  vector  along  the  drag  force  direction 

d  which  is  opposite  to  the  velocity  vector 

u. ,u  are  two  orthogonal  unit  vectors  defining  a  plane 

t  c  perpendicular  to  for  modeling  lifting  forces 

g  is  the  gravity  force  vector. 

A  state  model  including  3,  a  ,  and  a  as  state  variables  results 

U  G 

in  a  9-state  state  vector.  Further  discussion  of  this  problem 
is  found  in  Refs.  [7],  [8]. 

2.2.2  Filter  Gain  Computation 

The  filter  update  equation  (2.4)  gives  the  weighted 
sum  of  the  one-step  predicted  state  and  the  difference  of  the  new 
and  predicted  measurement  (this  difference  is  called  the  filter 
residual  process  and  in  more  rigorous  situations,  the  innovation 
process).  The  filter  gain  provides  the  weighting  of  this  update 
procedure.  When  the  filter  gain  is  very  small,  the  estimator 
becomes  insensitive  to  the  new  measurements  and  this  may  result 
in  large  bias  errors.  When  the  filter  gain  is  very  large,  the 
estimator  is  discarding  past  measurements  resulting  in  larger 
random  errors.  A  balance  can  be  achieved  if  one  has  sufficient 
knowledge  regarding  the  accuracy  of  the  state  model  and  the 
measurement  process.  Such  knowledge  is  embedded  in  the  filter  gain 
computation.  Four  approaches  for  treating  this  problem  are  discussed. 


2. 2. 2.1  The  Extended  Kalman  Filter  (EKF) 

If  an  extended  Kalman  filter  is  used,  the  filter 


gain  is 


T  -1 

K.  =  P,  ,.H,  R, 
k  k/k  k  k 


(2.12) 


where  is  the  measurement  Jacobian  matrix,  is  the  measure¬ 
ment  noise  covariance,  and  is  the  covariance  of  the  state 

estimate  obtained  by  solution  of  the  matrix  Riccati  equations. 


(1)  From  tk  to  tk+1,  compute  P^+jy^  by  integrating 


P  =  FP  +  PFa  +  Q  ,  P. 


(2.13a) 


(2)  At  tk+1 


Pk+l/k+l  Pk+l/k[I-Hk+l  (Hk+lPk+l/kHk+l  +Rk+1)  Hk+lPk+l/k} 

(2.13b) 

where  F  is  the  system  Jacobian  matrix  and  Q  is  the  process  noise 
covariance  matrix.  The  matrix  Q  must  be  chosen  to  be  representative 
of  modeling  errors.  For  example,  if  a  CV  model  is  chosen  while 
the  target  may  actually  be  accelerating,  Q  should  be  chosen 
proportional  to  the  expected  magnitude  of  target  acceleration. 

A  method  for  selecting  Q  such  that  the  filter  computed  covariance 
is  an  upper  bound  of  the  actual  performance  is  suggested  in  [9], 
Depending  upon  specific  applications,  the  value  of  Q  can  also 


be  obtained  empirically  using  Monte  Carlo  simulation.  Ref.  [8] 


Modeling  error  may  become  a  significant  contributor  to 
estimation  error  when  the  Riccati  equation  reaches  a  very  small 
steady  state  solution  (the  P  matrix) .  In  this  case  the  filter 
gain,  eq.  (2.12),  becomes  nearly  equal  to  zero  and  the  filter 
is  essentially  running  open-loop,  eq.  (2.4).  The  use  of  a  large 
Q  forces  the  steady  state  P  matrix  to  stay  "significantly" 
large,  so  that  the  estimate  is  sensitive  to  the  most  recent 
measurements . 

The  above  discussion  presents  a  heuristic  argument 
on  the  use  of  a  process  noise  covariance  Q.  Besides  methods 
using  covariance  upper  bounds  and  Monte  Carlo  studies  suggested  in 
[8]  and  [9],  there  also  exists  considerable  literature  which  dis¬ 
cusses  methods  of  estimating  Q,  R,  (and  sometimes  K  directly)  in 
realtime  based  upon  the  statistics  of  the  innovation  process.  One 
notable  paper  in  this  area  is  by  Mehra  [13],  We  briefly  outline 
Mehra's  method  below. 

It  is  well-known  that  if  both  the  system  model  and  noise 
statistics  are  true  representations  of  the  actual  physical  process 
the  filter  innovation  process 

-k  =  — k  ”  -(^k/k-l)  (2.14) 

is  white  Gaussian  with  zero  mean  and  covariance* 

♦This  is  clearly  not  true  for  nonlinear  systems.  In  most  applica¬ 
tions  using  an  extended  Kalman  filter,  this  expression  appears  to 
be  a  very  close  approximation. 


(2.15) 


HkP 


k/k-1 


\ 


and  the  filter  achieves  optimal  performance.  In  [13] ,  a  method 
for  testing  filter  optimality  based  upon  the  aforementioned 
property  of  the  innovation  process  was  presented.  If  the  test 
indicates  that  the  filter  does  not  attain  optimal  performance, 
one  then  proceeds  to  adjust  Q  and  R  so  that  the  covariance 
of  the  innovation  process  will  be  consistent  with  that  of 
filter  prediction.  This  method  first  computes  a  sampled 
correlation  function  assuming  that  yk  is  ergodic  over  a  certain 
time  interval,  one  therfore  has 


and 


Cj  £  E[Xi  Xi_jT]  =  HPwHT  +  R;  j=0 

=H[$>(I-KH)  ]j“14>[PooH'r-KC0]  ; 


(2.16) 
j  >o 


A 


N 


if  £  iA-j' 


(2.17) 


where  P^  denotes  the  steady  state  error  covariance  matrix  and  $ 
is  the  transition  matrix  of  the  system  dynamics. 

Using  the  above  equations  and  the  steady  state  Riccati 

equation,  Mehra  gives  a  procedure  for  solving  for  Q  and  R.  There 
are  situations  in  which  there  is  no  sufficient  number  of  indepen¬ 
dent  equations  for  solving  them?  Mehra  then  gives  a  recursive  pro¬ 
cedure  for  solving  for  the  filter  gain  K  directly. 
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We  emphasize  that  the  above  technique  is  important  be¬ 
cause  a  systematic  procedure  for  identifying  the  noise  covariance 

or  filter  gain  is  established.  It  may  not  be  very  useful  for 
realtime  applications  because  of  its  computational  requirement. 

For  non-realtime  applications,  however,  this  method  is  useful  for 
simulations  and  post-mission  data  analysis  studies. 

There  are  situations  in  which  the  dynamic  process  takes 
large  and  sudden  changes  (such  as  target  maneuvering) ,  such 
that  either  a  large  Q  must  be  used  all  the  time  to  account  for  the 
maximum  expected  deviation  or  else  a  quick  and  easy  method  must 
be  utilized  for  realtime  identification  of  Q.  This  subject  is 
covered  in  Section  3.1. 

There  exist  other  methods  for  preventing  the  P  matrix 
from  becoming  too  small.  These  including  the  finite  memory  filter 
[1] ,  [12],  and  the  fading  memory  filter  [1],  [10],  [11]. 

2. 2. 2. 2  The  Finite  Memory  Filter 

The  finite  memory  filter  applies  a  sliding  window 
to  the  data  and  computes  the  state  estimate  based  only  upon  data 
in  that  time  span.  The  window  width  is  selected  such  that  the 
system  model  is  an  adequate  approximation  to  the  actual  process 
over  the  time  interval.  During  this  time  interval,  the  system  is 
assumed  to  be  noise  free.  This  assumption  is  related  to  the 
selection  of  the  time  interval  in  that  the  variation  of  the 
unknown  parameter  over  this  interval  is  small.  The  filter  to  be 


discussed  is  due  to  Jazwinski  [12] .  We  briefly  outline  it  below. 
Let  the  measurement  sequence  be  denote  by 


—2 ' 


,-=m,'=m+l' 


'^k 


and  N=k-m  where  N  is  the  total  number  of  measurements  desired  in 
the  finite  memory  filter.  Let  the  finite  memory  state  estimate 
and  covariance  be  denoted  by  xk//m  k  and  P^ym  k,  respectively;  they 
can  be  computed  using  the  following  equations 

A  A  —  ]_A 

— k/m,k  =  Pk/m,k  (Pk/k  -k/k  “  Pk/m  -k/m) 


-1 


k/m,k 


=  P. 


k/k 


-1 


-  P. 


-1 


k/m 


(2.18) 


where  is  the  state  estimate  at  time  k  based  upon  all  data 

.  A  ,  , 

up  to  and  including  and  1S  the  state  estimate  at  time  k 

based  upon  all  data  up  to  and  including  z^.  One  can  interpret 
the  above  equations  to  say  that  the  finite  memory  estimate  is 
obtained  by  subtracting  an  estimate  based  on  all  data  prior  to 
the  time  window  from  the  estimate  based  upon  all  data.  We  note 
that  the  above  equations  can  be  re-organized  to  obtain  a  computa¬ 
tionally  more  efficient  and  numerically  more  stable  algorithm. 

For  more  discussions,  see  [12]. 

Notice  that  the  finite  memory  filter  requires  that  a 
batch  of  N  measurement  vectors  be  stored  (for  updating  . 

Furthermore,  the  computational  burden  is  much  larger  than  that  of 


13 


the  EKF.  Typically,  the  cost  of  added  computation  outweighs  the 
benefit  when  compared  with  other  methods. 

2. 2. 2. 3  The  Fading  Memory  Filter 

The  fading  memory  filter  (sometimes  referred  to  as  the 
aging  filter)  weights  recent  data  exponentially  higher  than 
past  data.  A  derivation  of  this  filter  can  be  found  in  [10]- [11]. 
The  resulting  algorithm  turns  out  to  be  very  simple.  Let  p£+jyk 
denote  the  covariance  of  the  one-step  predicted  estimate  of  the 
aging  filter,  it  is  related  to  the  Kalman  filter  covariance 
(solution  of  Eq.  (2.13a))  by 

Pk+l/k  =  aPk+l/k  (2.19) 

where  a  is  a  scalar  quantity  greater  than  unity.  The 
is  then  used  in  Eq.  (2.13b)  to  obtain  the  update  covariance 
Pk+l/k+l  The  scalar  a  -*-s  the  exponential  weighting  factor. 

This  is  accomplished  by  changing  the  measurement  noise  covariance 
matrix  to 

RJ  =  Ri(a)k"1  (2.20) 

for  i=l,...,k  where  k  is  the  current  time.  With  measurement 
noise  covariance,  data  rate  and  an  assumed  modeling  error,  one 
can  find  an  optimum  a  for  minimizing  the  mean  square  error. 
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2. 2. 2. 4  A  Constant  Gain  Filter:  The  u-B-y  Tracker 


In  some  cases,  because  of  computational  constraints, 
it  may  be  impractical  to  compute  the  filter  gain  in  real-time. 
Under  such  conditions  one  must  use  either  a  set  of  pre-computed 
filter  gains  or  a  constant  gain  filter.  One  commonly  used  con¬ 
stant  gain  filter  is  based  upon  the  steady  state  solution  of  the 
Riccati  equation.  In  order  to  minimize  the  modeling  error  effect, 
a  sufficiently  large  process  noise  covariance  Q  must  be  artfully 
selected. 

Another  popular  constant  gain  filter  is  the  oi-B-y 
filter,  (or  a- 1?  filter  when  using  a  CV  model),  [5],  [6].  The 

main  difference  between  the  a-B-y  filter  and  the  steady  state  gain 
filter  is  that  the  former  assumes  the  complete  independence  of 
the  three  spatial  coordinates  in  the  filter  update  equation. 

Notice  that  the  state  space  and  the  measurement  space  may  be 
related  through  a  nonlinear  function  (Eq.  (2.2)).  In  using  the 
a-0-y  filter,  this  assumption  is  not  allowed.  Let  (r,a,e) 
denote  the  radar  range,  azimuth  and  evaluation,  respectively; 
the  state  vector  becomes, 

,  •••  . ,T 

x  =  [r,r,r,a,a,a,e,e,e j  (2.21) 

Let  r  denote  [r,r,*r]T,  then  the  filter  update  equation  becomes 

— k+l/k+1  =  — k+l/k  +  Kk+1 k+1  “  -k+l/k]  (2.22) 

where  range  measurement  at  time  k+1.  The  gain  matrix 
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for  the  a-g-Y  filter  is 


K 


k+1 


Kk+1  *  fa'  t 


e 


2y 


(2.23) 


The  update  equations  for  azimuth  and  elevation  can  be  obtained 
accordingly. 

Notice  in  Eq.  (2.12),  when  an  extended  Kalman  filter  is 
used,  the  gain  matrix  has  the  dimension  (9x3)  while  the 
gain  matrix  for  the  a-g-Y  filter  consists  of  three  (3x1)  vectors. 

There  is  a  wide  range  of  methods  for  choosing  the 
values  of  a,g  ,  and  y.  One  method  is  to  compute  the  steady 
state  Kalman  gain  in  the  above  chosen  coordinate.  A  more 
reliable  method  is  to  conduct  extensive  Monte  Carlo  simulation 
studies  to  define  a,  g,  y  over  a  variety  of  cases  for  a  given 
application. 


2.3  A  Batch  Filter 

In  Section  2. 2. 2. 2,  we  discussed  a  finite  memory  filter 
developed  by  Jazwinski.  The  input  elements  of  that  algorithm  are 
the  outputs  of  the  recursive  Kalman  filter.  The  estimate  is 
based  upon  data  from  a  finite  time  interval  and  a  batch  of  N 
measurement  vectors  corresponding  to  that  time  interval  must 
be  stored.  Another  assumption  of  that  filter  is  that  the  system 
model  is  noise  free  and  this  assumption  is  reasonable  because 
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the  finite  time  interval  is  selected  to  reflect  the  tolerance  on 
model  errors.  The  Jazwinski  algorithm  is  optimum  for  linear 
systems  under  the  given  conditions. 

In  this  section,  we  present  an  algorithm  for  nonlinear 
systems.  This  algorithm  iteratively  solves  for  the  optimum 
estimate  based  upon  a  noise  free  system  model  and  a  batch  of  N 
measurement  vectors  z. z.2'***'^n*  filter  is  fundamentally 

different  from  the  approach  of  Jazwinski.  For  linear  systems, 
this  approach  is  the  same  as  the  Kalman  filter.  For  nonlinear 
systems  this  approach  gives  a  better  estimate  than  that  of  the 
extended  Kalman  filter. 

Let  the  system  and  measurement  equations  be  the  same 
as  eqs.  (2.1)  and  (2.2)  except  that  w  is  now  zero.  The  batch  of 
N  measurement  vectors  is  denoted  by  {z^,  z^,...,^}.  Let 
— N/N  denote  the  k“th  iteration  of  the  estimate  of  a^,  then 


Ak+1 

-N/N 


^N/N 


+  P 


N/N 


r  N 

to* 

L^i  n 


’  kT  -l 
H*  R  1 (z 
n  n  — n 


(2.24) 


[EgJVVwJ 
LnTl  n  n  n  n  nJ 


(2.25) 


where 


n=N- 1 ,  N-2 , . . . ,  1 
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-n/N 


I  (an  identity  matrix) 


the  Jacobian  matrix  of  f (xg^N) 


the  Jacobian  matrix  of  h(x^N) 


result  of  integrating  x=f(x)  backward  from 
tN  to  tn  using  ^  =  x^/N. 


''k-f-l  Ak 

The  iteration  is  terminated  when  and 


are  sufficiently 


close. 


We  make  the  following  remarks: 

The  above  algorithm  is  a  realization  of  the  maximum 
likelihood  estimator  with  Gaussian  measurement  noise  process. 
It  is  well-known  that  the  maximum  likelihood  estimate 
is  asymptotically  efficient  and  Gaussian  and  approaches 
the  Cramer- Rao  bound. 


The  PN/N  of  (2.25)  is  an  approximate  expression  for 


the  covariance  of  xN/N» 


The  PN/N  evaluated  at  the  true 


state  is  the  Cramer-Rao  lower  bound  on  the  covariance 

A  A 

of  Since  approaches  the  true  state 


since  Sn/n 


with  probability  one,  PN/N  also  approaches  the 
Cramer-Rao  bound  with  probability  one. 


Notice  that  the  inverse  of  P. 


is  Fisher’s  informa¬ 


tion  matrix.  The  invertibility  of  the  information 

matrix  is  tied  with  the  observability  of  the  system, 
see  for  example  [58] . 


4)  For  linear  systems  a  closed  form  solution  can  be  found 
and  the  iterative  procedure  becomes  unnecessary.  An 
interesting  exercise  is  to  derive  a  batch  filter  using 
the  polynominal  dynamics  with  linear  measurements  equa¬ 
tions  (see  for  example  [49] ) . 


There  are  many  application  areas  for  this  algorithm. 
For  example,  in  tracking  space  objects  where  the  target  dynamics 
can  be  modeled  very  accurately,  the  algorithm  of  this  section  is 
particularly  suitable.  This  method  has  been  used  for  ballistic 
trajectory  tracking  with  angle-only  measurements  [30]  and 
tracking  of  deep  space  satellites,  [31],  Another  application 
is  for  track  initiation .  Since  the  initial  covariance  and  state 
estimates  are  not  generally  given  a  priori,  the  above  algorithm 
can  obtain  the  best  estimates  based  on  the  first  N  measurement 

/V 


vectors  and  then  proceed  to  use  x^/N  anc*  PN/N  as  t^ie 

state  and  covariance  estimates,  respectively.  This  method  is 

sometimes  referred  to  as  the  information  matrix  approach  for 


filter  initiation. 


2.4  Summary 

In  the  above,  we  have  discussed  various  algorithms  for 
addressing  the  basic  tracking  problem.  These  approaches  employ 
simple  to  sophisticated  system  models  for  target  dynamics  and 
attempt  to  compensate  for  modeling  errors  in  a  variety  of  ways. 
Several  well-known  nonlinear  estimation  algorithms  were  not  dis¬ 
cussed.  These  include  the  second  order  filter  [14]  and  the  single 
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stage  iterative  filter  [15] .  They  are  not  included  because  of 
their  excessive  computational  requirement  although  these  algorithms 
can  indeed  improve  the  estimation  accuracy. 


Before  closing  this  section,  a  brief  algorithm  comparison 
can  be  stated.  If  the  objective  of  tracking  is  to  obtain  pre¬ 
cision  information  about  the  target  dynamics ,  then  one  should  use 
the  most  accurate  target  model  and  apply  the  EKF,  (or  even  more 
sophisticated  algorithms).  If  the  dynamic  model  is  sufficiently 
accurate  so  that  the  process  noise  term  is  negligible,  then  the 
algorithm  considered  in  the  Section  2.3  is  a  good  choice  provided 
that  the  computation  time  and  data  storage  requirements  are  not 
excessive. 

If  the  objective  is  just  to  maintain  the  target  in  track, 
then  one  may  use  the  simplest  track  algorithm  such  as  the  a-B-y 
tracker.  One  exception  to  this  case  is  when  tracking  in  a  dense 
target  environment  where  precision  tracking  may  be  necessary  for 
target  correlation.  This  subject  will  be  discussed  later. 

The  finite  and  fading  memory  filters  are  usually  the 
secondary  choices  (especially  the  finite  memory  filter)  because 
the  same  purpose  (reducing  sensitivity  to  model  errors)  can  be 
achieved  by  adjustment  of  the  process  noise  covariance,  Q.  They 
are  nevertheless  included  here  because  (1)  they  reflect  the  histori¬ 
cal  development  of  adaptive  filtering  techniques  and  (2)  they  still 
provide  an  option  for  readers  to  choose  for  their  applications. 
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3. 


TARGETS  WITH  SUDDEN  MANEUVERS 


Targets  with  sudden  maneuvers  can  be  modeled  as  systems 
with  abrupt  changes.  We  modify  equation  (2.1)  to  become  two  sets 
of  equations,  one  representing  the  pre-maneuver  dynamics  and  the 
other  incorporating  the  maneuver  feature 

x  =  f  (x,jO  +  w  (2.1a) 


where  x  is  the  vector  representing  maneuvering  force  and  satisfies 
m 


X  =  0 

— m  — 


for  t  >  t 


(3.1) 


=  Ln  +  f°r  t  >  t_ 

— m  — m  — m  — m  —  m 


(3.1a) 


where  tm  is  the  time  the  maneuver  begins,  f^t* )  is  the  maneuvering 
dynamics  and  is  the  system  noise  for  )  •  For  targets  with 

sudden  maneuvers,  tm  is  unknown,  (• )  and  may  be  unknown  or 
partially  known. 

In  tracking  airplanes  for  example,  the  target  may  first 
be  flying  in  a  straight  line  with  constant  speed.  A  CV  model 
is  adequate  for  the  airplane  dynamics  in  this  case.  A  sudden 
maneuver  of  this  airplane  implies  that  the  airplane  is  accelera¬ 
ting  unexpectedly  and  the  acceleration  is  time-varying  and  follow¬ 
ing  an  unknown  profile.  The  instantaneous  acceleration  vector  is 
therefore  the  maneuvering  vector  x^  defined  above. 

In  tracking  a  ballistic  re-entry  vehicle,  the  target 
dynamics  is  the  equation  (2.11)  with  a  and  a  equal  to  zero.  A 

C  t 
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sudden  aerodynamic  maneuver  means  ac  and  at  become  nonzero 

and  follow  unknown  time-varying  profiles.  The  maneuvering 

vector  in  this  case  contains  a  and  a .  as  its  elements. 

c  t 

Four  approaches  to  this  problem  are  discussed  individually 

below. 

3.1  Filter  Compensation  Using  Process  Noise  Covariance 

In  this  first  method,  one  simply  ignores  the  maneuver 

vector  x  and  lunps  the  system  errors  introduced  by  x  with  the 
— m  — m 

process  noise  term  w^.  If  the  estimator's  only  concern  is  to 
maintain  the  target  in  track  (adequate  position  estimation 
accuracy) ,  this  method  can  work  quite  well. 

Basically,  it  examines  the  •'regularity"  of  the  filter 
residual  vector 

Xk  -  £k  -  H<£k/k-i>  (3.2) 

against  its  covariance  matrix 

\  -  HkPk/k-lHkT  +  \  (3*3) 

using  (the  Chi-square  variable) 

**  *  *\\  ■  ,3’4) 

When  Jtk  becomes  too  large  one  suspects  that  the  target  is  man¬ 
euvering  and  the  covariance  of  w,  Q  is  increased  so  that  is  re¬ 
duced  to  a  reasonable  value.  This  method  therefore  has  the 


22 


combined  feature  of  maneuver  detection  and  filter  compensation. 

This  method  is  based  upon  the  adaptive  filter  of 
Jazwinski,  [16]  and  shares  great  similarities  with  the  method 
of  Mehra,  [13]  (see  also  section  2. 2. 2.1).  The  main  difference 
is  that  the  above  method  tests  the  filter  regularity  using 
Chi-square  statistics,  and  stresses  the  simplicity  (not  the 

optimality)  of  the  approach. 

A  thorough  discussion  of  this  method  and  its  performance 
against  maneuvering  re-entry  vehicles  can  be  found  in  Refs.  [7] 

[8], 


3.2  State  Augmentation 


The  second  method  is  straightforward,  computationally 
more  costly,  but  with  substantially  better  performance  than  the 
previous  method.  This  method  is  to  include  as  part  of  the 
state  vector,  i.e.,  the  augmented  state  consists  of 


2a  =  ^  O'* 


(3.5) 


In  the  case  when  the  target  maneuvering  dynamic  is  completely 
unknown,  one  uses 


2m 


=  w_ 

“HI 


(3.6) 


where  w^  is  modeled  as  a  zero  mean  white  noise  process  with 
covariance  Qm.  If  bounds  on  the  magnitudes  of  maneuvers  are 
known,  a  method  for  choosing  Qm  such  that  the  actual  filter  per- 
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formance  is  bounded  by  the  computed  filter  covariance  can  be  found 

in  [91.  As  a  rule  of  thumb,  the  values  of  the  entries  in  Q  should 

m 

be  a  fraction  of  the  expected  magnitude  squared  of  the  maneuver 
force  and  proportional  to  the  measurement  time  interval.  A  method 
for  selecting  Qm  for  estimating  ballistic  and  lifting  coefficients 
of  reentry  vehicles  is  discussed  in  [7],  f 8 ] . 

Notice  that  Eq.  (3.6)  assumes  uncorrelated  in  time. 
This  assumption  is  sometimes  not  very  realistic.  A  model  often 
used  in  airplane  tracking  is 


=  (txm 
— m  — m 


a  =  w 


(3.7) 


where  a  is  the  correlation  constant  to  be  estimated  and  is 
a  noise  process.  Methods  for  selecting  statistics  for  w^  and  the 
performance  against  airplane  tracking  can  be  found  in  Ref.  [17J. 

The  maneuvering  state  x^  is  usually  not  influenced  by 
the  state  vector  x  (Eqs.  (2.1a)  and  (3.1a)).  With  this 
assumption,  one  can  compute  the  state  and  maneuver  estimates 
separately  to  obtain 


-k/k  “  -k/k  +  \  2m, k  (3*8) 

where  x  is  the  estimate  assuming  x  is  zero,  x  ,  is  the  maneuver 
—  — m  — m,K 

estimate,  is  a  gain  matrix,  and  x^^  is  the  final  state 
estimate.  The  advantage  of  this  decoupled  implementation  is  a 
saving  in  computation.  It  can  be  shown  that  the  decoupled  estimator 


is  optimum  for  linear  systems  when  tm  is  known  and  the  maneuver 
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4 


state  model  is  known  and  deterministic;  see  for  example  [18],  [19], 
and  [21].  Suboptimal  designs  for  the  case  when  the  above  assump¬ 
tions  are  not  true  can  be  found  in  [9] . 

Since  the  maneuver  time  t  is  unknown,  a  trival 

m 

application  of  the  above  method  is  to  use  it  throughout  the  entire 
track.  Inevitably,  the  filter  performance  is  degraded  when  the 
target  is  non-maneuvering,  Refs.  [8],  [9],  An  ideal  extension  is 
to  use  a  maneuver  detector  for  switching  the  tracking  filter  from 
nonmaneuver  to  maneuver  mode,  this  subject  is  discussed  in  the 
following  subsection. 

3 . 3  Manuever  Detection 

The  detection  problem  is  the  problem  of  discriminating 


the  following  two  hypotheses  based  upon  filter  residuals: 


H^s  Xk  =  Xk  +  9^ (xm) '  maneuvering  target  hypothesis 


H0!  *k  =  Ik 


;  non-maneuvering  hypothesis 


(3.9) 


for  k  =  kg,...,  K  where  is  the  residual  vector  before  the  maneuver 
starts,  and  g(  )  is  a  known  function  relating  the  maneuvering  vector 
3^  to  the  residual  vector.  A  generalized  likelihood  ratio  test 


max  ~  . 

A  =  2m  P(Xk;k0,...>K/H1)  >  ^ 

P  <Xk' ^0' * • * ' R/Ho) 


(3.10) 


For  a  given  application  (target  dynamics,  sensor  type,  etc.),  the 
above  equation  can  be  further  simplified,  [7],  [18],  [20],  [21]. 

The  use  of  a  generalized  likelihood  ratio  test  for 
detecting  maneuvers  is  discussed  in  [7],  [8]  for  re-entry  vehicles 
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and  in  [20]  for  airplanes.  General  discussions  on  the  detection  of 
sudden  changes  in  linear  systems  can  be  found  in  Refs.  [18]  and 
[21]  . 


Drawbacks  of  using  a  detector-directed  tracking  include 
(1)  detection  delay  and  (2)  large  transient  errors  during  filter 
switch.  Numerical  results  in  Ref.  [8]  show  large  estimation 
errors  immediately  after  filter  switch.  This  is  partially  due 
to  the  maneuvering  filter  going  through  its  transient  period. 

Notice  that  the  maneuver  detection  relys  on  the  fact 
that  target  maneuvering  generates  residual  bias  in  a  non-maneuver¬ 
ing  tracker.  Once  the  filter  is  switched  to  maneuvering,  the 
above  detection  scheme  can  not  be  used  to  discriminate  if  the 
target  has  returned  to  non-maneuvering  status. 

All  the  above  problems  can  be  alleviated  if  one  employs 
the  adaptive  multiple  model  estimator,  [22]  —  [29],  discussed  in  the 
next  subsection. 


3.4  Multiple  Model  Estimator 


Let  denote  the  hypothesis  that  the  target  is  maneuver¬ 
ing  and  Hq  the  hypothesis  that  the  target  is  non -maneuvering. 

One  may  construct  two  filters,  the  first  one  uses  a  maneuvering 


dynamic  while  the  other  one  uses  a  non-maneuvering  model.  Let 
ve  state  be  denote 
is  obtained  using 


these  respective  state  be  denoted  by  and  x^k;  the  optimum 


estimate 


-k/k  *  Pl(k)£k/k  +  P0(k)-k/k 


(3.11) 


where  (k)  is  the  a  posteriori  probability  that  is  true  at 
time  t^.  The  computation  of  Pi(k)  can  be  very  complicated.  If 
one  uses  a  simplifying  assumption  that  hypotheses  at  time  tk  are 
independent  of  hypotheses  before  t^  i.e.,  a  first  order  Markov 
process,  then  one  obtains  a  much  simplified  expression  for  P^(k). 


Epij(k)pijpj(k_1) 


P^k)  = 


(3.22) 


E  E  P..(k)p.  p  (k-D 

1=0  j=0  J  J 

where  p^^ (k)  is  the  residual  density  assuming  that  the  i-th 

hypothesis  is  true  at  time  t^  and  the  j-th  hypothesis  is  true  at 

time  t,  1 ,  P.  .  is  the  transition  probability,  and  P .  (k-1)  is  the 
k-i  j 

a  posteriori  probability  that  H.  is  true  at  time  t^^.  Notice 
that  we  have  made  the  above  equation  slightly  more  general  by 
assuming-that  there  are  total  of  N+l  hypotheses.  The  xj^k  is 

computed  using 


4/k 


(3.23) 


P(H.(k-l)/Hi(k),Zk)  = 


Pij  (k)pijpj (k-1) 


E  p. 


(3.24) 
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where  x^y^  is  the  state  estimate  using  updated 

with  the  i-th  hypothesis  and  denotes  all  the  measurements  up  ^  time  t^. 

Derivations  and  discussion  of  the  above  results  in  a 
more  general  context  are  given  in  [22]  ,  [23] . 

The  advantage  of  the  above  approach  is  that  when  the 
target  switches  between  maneuvering  and  non-maneuvering  modes,  the 
hypothesis  probability  values  change  to  provide  a  smooth 

A 

transition  of  the  final  estimates  x^^. 

A  higher  level  multiple  model  approach  is  to  use  several 
hypotheses  to  model  different  maneuvering  force  levels.  This 
approach  enhances  the  estimation  performance  at  the  cost  of  more 
computation  resources. 

We  would  like  to  emphasize  that  the  multiple  model 
estimator  is  a  general  adaptive  estimation  technique.  It  was 
first  derived  by  Magill  [27]  for  the  time  invariant  hypothesis 
case,  i.e., 

/ 

1  ,  when  i  =  j 


P.  .  = 
13 


(3.25) 


otherwise 


Its  extension  to  the  switching  hypothesis  case  was  the  subject 
of  [22].  References  [28]  and  [29]  stated  the  above  concept  in  the 
continuous  domain  and  gave  a  representation  theorem  known  as  the 
partition  theorem. 

A  word  of  caution  about  the  multiple  model  method;  it 
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is  a  tightly  tuned  algorithm  (for  the  case  of  modeling  various 
maneuvering  levels)  and  therefore  very  vulnerable  to  the 
situation  where  none  of  the  models  match  with  the  actual  dynamics. 
In  this  case,  strange  filter  behavior  may  appear  and  this  behavior 
is  not  yet  completely  understood  from  a  theoretical  point  of  view, 
[23],  [24]. 

Discussions  applying  this  method  to  target  tracking 
can  be  found  in  [25]— [26] .  Its  extension  to  adaptive  control  was 
applied  to  the  flight  control  of  a  F-8C  experimental  aircraft  [24], 
A  tutorial  treatment  of  this  method  for  state  estimation  is  given 
in  [23], 


3.5  Summary 

In  this  section,  we  have  discussed  four  approaches  to 
tracking  targets  with  sudden  maneuvers.  The  first  method 
(Section  3.1)  provides  a  way  for  adjusting  the  process  noise 
covariance  level  through  use  of  the  filter  residual  "regularity" 
to  compensate  for  the  modeling  error  induced  by  target  maneuver. 
This  method  uses  the  least  computation  but  does  not  give  very 
precise  velocity  and  parameter  estimates.  The  second  method 
(Section  3.2)  is  to  augment  the  state  vector  with  maneuvering 
variables.  The  dimension  of  the  state  vector  is  enlarged  and  the 
filter  is  therefore  computationally  more  costly.  It  does  however 
provide  more  accurate  state  estimates.  The  drawback  to  this 
method  besides  the  higher  computational  burden  is  that  the  filter 
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performance  is  degraded  when  the  target  is  non-maneuvering.  A 
method  to  circumvent  this  problem  is  to  use  a  maneuver-detector- 
directed  filter  (Section  3.3).  This  method  suffers  from  large 
errors  occurring  during  filter  switch  over.  A  compromise  to  this 
problem  is  to  contruct  two  filters  (Section  3.4)  with  one  using 
a  maneuvering  dynamic  model  and  the  other  a  non-maneuvering 
dynamic  model.  The  final  estimate  is  a  weighted  sum  of  outputs 
of  these  two  filters  using  the  a  posteriori  hypothesis  probabili¬ 
ties  as  weighting  factors.  This  method  has  the  advantages  of  all 
the  above  approaches  but  at  a  cost  of  a  much  larger  computational 
burden.  One  may  also  use  a  bank  of  filters  to  model  different 
maneuver  levels  and  apply  the  multiple  model  adaptive  estimation 
method.  This  method  may  work  very  well,  it  is  however,  also  very 
sensitive  to  model  mismatch  errors  and  system  nonlinear ites. 


4. 


TRACKING  WITH  ANGLE-ONLY  MEASUREMENTS 


In  this  section,  the  problem  of  tracking  with  a  sensor 
measuring  only  the  target  line-of-sight  angle  is  discussed. 

There  may  be  two  separate  objectives  of  tracking.  The 
first  one  is  simply  trying  to  maintain  targets  in  track  in  the 
angle  domain.  In  this  case,  techniques  discussed  in  Section  2 
are  applicable.  The  only  difference  is  that  ir  the  angle  only  case 
the  target  dynamic  equations  are  described  in  two-dimensional 
coordinates.  Polynomial  equations  decoupled  in  two  orthogonal 

angular  directions  are  often  used.  The  second  tracking  objective 
is  to  obtain  estimates  of  the  complete  state  vector  as  defined  in 
a  three-dimensional  coordinate  system.  This  objective  may  not 
always  be  achievable  since  tracking  with  angle-only  measurements 
may  constitute  an  unobservable  system.  Physically,  it  can  be 
explained  as  follows.  In  radar  tracking  systems,  each  measurement 
vector  determine  the  instantaneous  target  position  to  within  a 
finite  uncertainty  volume,  i.e.,  uncertainties  in  both  range  and 
angles  can  be  expressed  with  finite  standard  deviations.  In  an 
angle-only  tracking  system,  the  uncertainty  volume  of  each  measure¬ 
ment  vector  is  infinite  (due  to  the  inability  to  measure  range) . 
Such  a  system  may  be  observable  only  for  certain  types  of 
target  dynamics.  For  example,  when  a  telescope  is  used  to  track 
a  satellite,  this  constitutes  an  observable  system  and  target 
range  can  be  estimated  (with  large  errors  however)  because  the 
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satellite  trajectory  is  influenced  by  the  earth  gravity  [30], 

If  a  telescope  is  used  to  track  an  airplane  traveling  with  constant 
speed,  this  system  is  unobservable,  i.e.,  the  estimation  of  the 
3-dimensional  dynamics  is  impossible. 

There  are  means  available  for  enhancing  the  observabil¬ 
ity  of  the  system.  For  example,  one  may  use  two  passive  sensors 
at  separate  locations  simultaneously  tracking  the  same  object. 

The  intersection  of  two  angular  beams  gives  the  total  measurement 
uncertainty  which  now  has  a  finite  uncertainty  volume.  One 
complication  of  this  approach  occurs  when  tracking  in  a  multiple 
target  environment.  Here,  one  is  confronted  with  the  problem  of 
recognizing  the  same  target  at  both  sensors.  The  subject  of 
sensor-to-sensor  correlation  and  efficient  methods  of  processing 
multiple  sensor  measurements  are  discussed  in  the  next  section. 
Another  method  of  enhancing  observability  is  to  incorporate  other 
types  of  measurements  in  addition  to  angle  measurements.  One 
such  applications  is  in  passive  sonar  tracking  [32]  where  target 
Doppler  measurement  is  included.  A  passive  sonar  system  with 
Doppler  measurements  makes  the  system  completely  observable. 

We  note  that  orbit  estimation  using  angle  only  measure¬ 
ments  is  an  ancient  problem  (see  [33]  for  a  discussion  of  classical 
approaches  and  a  list  of  references).  Most  emphasis  of  [33]  how¬ 
ever,  was  placed  on  the  plant  tary  mechanics  and  dynamic  modeling 
and  very  little  effort  was  applied  from  the  estimation  point  of  view. 
For  example,  the  Gauss  method  for  orbit  determination  ensures  the 


solution  of  the  entire  orbit  with  only  three  angle  measurements. 
This  is  true  if  noise  free  angle  measurements  can  be  made.  With 
even  slight  errors  in  angle  measurements,  the  target  range 
estimation  error  can  be  very  large. 

When  one  is  tracking  a  ballistic  object  (a  satellite 
or  a  long  range  ballistic  missile)  using  a  passive  sensor,  one 
is  confronted  with  two  questions:  1)  how  to  initiate  a  Kalman 
type  recursive  tracking  filter;  and  2)  since  the  range  estimation 
error  will  inevitably  be  larger,  will  the  extended  Kalman  filter 
provide  adequate  performance.  In  [30],  these  issues  were  studied 
in  detail.  It  was  shown  that  the  batch  filter  described  in 
Section  2.3  can  be  used  to  provide  initial  conditions  for  the 
extended  Kalman  filter.  Other  discussions  include  methods  for 
computing  an  initial  guess  for  the  iterative  procedure,  applica¬ 
tions  of  the  batch  filter  recursively  for  tracking,  and  techniques 
for  incorporating  trajectory  a  priori  knowledge  (bounds  on  heading 
angles,  velocity,  energy,  etc.)  for  improving  the  state  estimation 
accuracy.  An  important  conclusion  of  [30]  is  that  the  performance 
of  the  batch  filter  asymptotically  approaches  the  Cramer-Rao 
bound  (see  Section  7.3  of  this  paper)  for  the  covariance  of  state 
estimates  while  the  performance  of  the  extended  Kalman  filter 
generally  does  not  have  this  property. 

Finally,  we  remark  that  recent  developments  in  target 
tracking  on  the  focal  plane  using  a  infrared  sensor  are  in  the  pro- 
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blem  areas  of  both  signal  processing  and  target  tracking.  The 
approach  discussed  in  [55]  is  to  first  process  focal  plane  data 
with  a  prescribed  model,  then  establish  two-dimensional  target 
state  estimates  using  a  Kalman  filter. 


i 


m 
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5.  TRACKING  WITH  MULTIPLE  SENSORS 

There  are  a  number  of  applications  in  which  multiple 
sensors  operate  in  a  multiple  target  environment.  There  are 
two  problems  which  confront  a  multiple  sensor  tracking  system: 

1)  how  to  efficiently  update  the  tracking  filter  with  multiple 
sensor  measurements  and  2)  how  to  identify  the  same  target  for 
all  sensors  in  a  multiple  target  environment.  The  first 
question  is  discussed  in  Section  5.1.  The  second  question  is 
a  special  case  of  tracking  in  a  dense  target  environment.  It 
fits  however,  conveniently  in  the  context  of  sensor-to-sensor 
correlation  and  will  therefore  be  discussed  in  Section  5.2. 

5.1  Filter  Update  Algorithm  Considerations 

Consider  the  situation  in  which  there  are  several  sen¬ 
sors  simultaneously  observing  a  single  target.  If  these  measure¬ 
ments  are  not  time  synchronized,  one  can  formulate  the  problem 
as  that  of  state  estimation  with  non-uniform  measurement  times. 

It  should  be  noted  that  in  this  particular  formulation  the 
functional  transformation  from  state  space  to  observation  space 
depends  upon  the  geometry,  hence  one  can  be  confronted  with  a 
situation  in  which  the  measurement  equations  change  from  observa¬ 
tion  to  observation. 
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If  the  measurements  are  time  synchronized  (such  as  in 
a  multi-static  radar  system) ,  then  one  is  interested  in  seeking 
the  most  efficient  way  of  processing  these  measurements.  Let  z^  ^ 
denote  the  measurement  at  time  t^  from  the  i-th  sensor,  then 

zk  ^  “  h^  (x^ )  +  vk  l  *  I#  (5.1) 


Three  processing  methods  are  discussed  by  Willner,  et  al.,  [34], 


1.  Parallel  filter:  The  set  of  measurement  vectors  are 
concatenated  to  form  a  new  measurement  vector 


2. 


r2T  -T  iT 

— k ,  1 '  — k ,  2 '  *  *  * '  — k ,  I J 


(5.2) 


which  is  processed  by  the  filter  at  once. 

Sequential  filter:  Each  measurement  is  processed 
sequentially  by  the  filter  with  zero  prediction  time 
between  measurements. 


3.  Data  compression  filter:  Prior  to  processing, 
measurement  vectors  are  transformed  to  a  common 
coordinate  system  and  compressed  to  a  single  measure¬ 
ment  vector.  Let  z>  i  denote  the  measurement  of  i-th 
sensor  and  all  have  been  transformed  to  a 

common  coordinate  system.  The  compressed  measurement 
vector  z,  is 


^k  =  *k 


i=l 


,i  -k,i 


(5.3) 


(5.4) 
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The  above  three  methods  can  be  shown  to  be  equivalent 


and  optimum  for  linear  systems,  [34],  The  data  compression  method 
is  computationally  the  most  efficient  for  a  wide  range  of  cases  [34] . 

It  is  possible  to  use  decentralized  processing  wherein 
each  sensor  produces  its  own  estimate.  These  are  then  combined  j 

to  form  a  single  estimate.  The  major  drawback  to  such  an  approach 
is  that  the  computational  burden  is  significantly  higher  than  any 
of  the  methods  described  above.  Furthermore,  when  the  system  model 
is  stochastic,  the  state  estimates  obtained  by  each  sensor  are 
correlated.  It  is  uncertain  as  to  whether  the  estimate  combination 
method  is  optimum  even  if  the  correlation  matrics  in  the  combina¬ 
tion  procedure  are  included  (see  for  example,  the  Gauss-Markov 
theorem  stated  in  [35]). 

Algorithms  for  data  compression  with  non-synchronized 
data  vectors  is  a  straightforward  extension  of  the  work  described 
above.  This  may  involve  integration  of  measurement  vectors  to 
achieve  time  alignment  followed  by  an  application  of  the  data 
compression  equations.  Polynomial  data  smoothing  may  be  used  for 
data  compression  in  very  high  data  rate  systems  as  is  discussed  in 
[36],  One  can  also  apply  polynomial  smoothers  to  the  non- 
synchronized  data  case. 

5.2  Sensor-to-Sensor  Correlation 

One  important  problem  facing  a  multiple  sensor  tracking  i 
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system  in  a  multiple  target  environment  is  the  unique  identifica¬ 
tion  of  the  same  target  as  observed  by  each  sensor.  There  are 
two  approaches  to  this  problem.  The  first  approach  attempts  to 
correlate  the  existing  track  files  (state  estimate  of  a  given 
string  of  measurements)  with  measurements.  This  method  is  the 
same  as  the  track  continuation  mode  of  tracking  in  the  dense 
target  environment.  We  will  address  this  approach  in  the  next 
section.  The  second  method  is  to  directly  correlate  the  set  of 
measurements  from  the  i-th  sensor  with  that  of  the  j-th  sensor. 

A  brief  discussion  of  the  second  method  is  given  below. 

Consider  the  case  when  there  are  two  sensor  simultane¬ 
ously  tracking  multiple  targets.  The  results  which  follow  can 
be  extended  to  that  of  multiple  sensors.  Let  i=l,...,N} 

and  {z^;  j=l,...,N}  denote  N  measurements  obtained  by  the  first  and 

* 

second  sensors,  respectively  and  respective  covariances  are  given 
by  and  E.. .  Consider  now  the  question  of  which  ^  corresponds 
to  a  particular  z^.  This  is  an  N-array  hypothesis  decision  pro¬ 
blem.  That  is,  for  a  given  z_. one  asks  which  of  the  following 
hypotheses  is  true. 


*We  have  obviously  neglected  the  problem  of  imperfect  detection  at 
the  sensor  and  unequal  coverage  problem  by  assuming  that  each  sensor 
observes  the  same  number  of  targets  and  that  there  is  a  one-to-one 
assignment  between  them.  For  dicussions  which  include  these  fac¬ 
tors,  consult  Ref.  [37], 
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where  n.  is  a  random  noise  vector  with  zero  mean  and  covariance 

—l 

E^.  Using  the  likelihood  ratio  test  procedure,  it  can  be  shown  that 
the  is  choosen  as  the  one  which  maximizes  the  joint  density  of 
z ^  and  y. j  •  That  is,  decide  that  is  true  when 

=  max  (5.6) 


If  z.  and  y-  are  Gaussian  random  vectors,  then  one  obtains  the 


-i  *0 

following  equivalent  procedure. 


Decide  is  true  when  Yj  minimizes 
"if  “  <ii-i.,,T(pi+E«>':L(2i-^) 


(5.7) 


for  all  £=1,...,N.  Notice  that  (5.7)  is  a  Chi-square  random 
variable  or,  a  weighted  distance  measure. 

The  above  discussion  gives  a  procedure  for  selecting 
a  Yj  vector  for  a  given  z^  vector.  When  one  considers  all  z^ 
vectors,  this  procedure  can  not  be  extended  without  modif icantions. 
This  is  because  if  one  repeats  the  above  procedure  for  all  z^’s,  one 
may  obtain  the  correlation  of  the  same  measurement  from  one  sensor 


to  several  measurements  of  the  other  sensor.  One  therefore  has 

to  impose  the  constraint  that  each  measurement  can  only  be  assigned 

(correlated)  once  while  optimizing  some  performance  index  (the 

problem  of  multiple  correlation  due  to  limited  sensor  resolution 

is  discussed  in  [37]).  One  such  performance  index  is  the  sum  of 
* 

all  Chi-squares.  That  is,  the  N  correlated  pairs  of  and  ^ 
are  those  achieving  the  minimum  of 


under  the  constraint  that  each  z ^  and  ^  can  only  be  used  once. 

We  note  that  the  above  problem  is  the  same  as  the 
assignment  problem  in  operations  research.  The  optimum  answer  may 
be  obtained  by  exhaustively  searching  for  all  combinations  which 

results  in  searching  through  N!  possibilities.  A  procedure 
called  the  Hungarian  method  [33]  (or  the  Munkre's  method  for  a 

specific  implementation  procedure  [39])  requiring  at  most 
3  2 

(11  N  +  12  N  +  31  N)/6  operations  is  often  used  in  this  type  of  problems. 

In  the  case  when  the  target  density  is  not  very  high  so 
that  the  w^j  of  a  mismatched  pair  may  attain  large  values,  a  thres¬ 
hold  may  be  first  applied  to  examine  all  w^'s.  For  Gaussian 
measurement  vectors,  this  threshold  can  be  selected  to  provide  a 
given  probability  of  leakage.  Those  pairs  exceeding  the  thres¬ 
hold  are  first  rejected,  one  is  therefore  only  required  to  corre¬ 
late  on  a  subset  of  measurement  pairs  which  do  not  exceed  the  threshold. 

♦This  criterion  is  corresponding  to  the  minimum  error  probability 
for  Gaussian  measurement  errors. 
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6.  TRACKING  IN  A  DENSE  TARGET  ENVIRONMENT 


Tracking  in  a  dense  target  environment  was  the  subject 
of  a  recent  review  article  [40]  and  several  invited  sessions  in 
recent  IEEE-CDC  conferences.  Representative  work  in  this  area 
can  be  found  in  [40]- [48].  This  problem  is  sometimes  referred  to 
as  scan-to-scan  correlation  or  tracking  data  association.  In 
this  section,  we  offer  some  very  general  discussion  on  this  pro¬ 
blem  and  the  interested  reader  should  consult  the  references  for 
details. 

This  problem  can  be  divided  into  two  phases.  The  first 
phase  is  track  initiation  and  the  second  phase  is  track  main¬ 
tenance.  They  are  discussed  individually  below. 

6.1  Track  Initiation 

Consider  the  case  of  a  scanning  sensor.  The  first  and 
second  scan  produce  N^  and  N2  detections,  respectively.  The 
problem  is  to  associate  the  two  sets  of  detections  to  form  min 
(NlfN2)  number  of  track  files.  Notice  that  we  have  assumed  that 
N1  ^  N2*  Th*s  can  be  c®use^  by  (1)  imperfect  detection  (2) 
emergence  of  new  targets  in  the  second  scan,  and  (3)  targets 
leaving  the  sensor  field  of  view  before  the  second  scan.  In  the 
following,  an  approach  for  track  initiation  with  k  scans  of  data 
is  described. 

Let  denote  all  the  measurements  (N)  collected  during 
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k-th  scan,  i.e., 


Z^  —  { 2^  ( k) ,  z 2  ( k)  ,  • . •  ,  2_jj  ( k)  } 


(6.1) 


Let  Z  denote  the  set  of  measurements  up  to  and  including  the  k-th 


scan,  i.e. 


{ Z± ;  i— 1,...,  k} 


(6.2) 


Furthermore  to  simplify  discussion,  assume  that  N  is  the  number  of 
detections  for  all  Z  '  s.  Assume  also  that  sensors  have  perfect  tar¬ 
get  detection.  When  this  is  not  true,  one  has  to  enumerate  more 

1c 

hypotheses  to  account  for  all  possibilities.  With  Z  ,  there  can  be 
N  combinations  of  measurement  sequences  and  each  measurement  se¬ 
quence  represents  a  possible  track.  Let  each  possible  combination 
be  denoted  by  a  hypthesis,  Hm  (k)  which  is  defined  by 

K 


H  (k)  =  {z  (1),  z  (2) , ... ,z  (k) > 
mk  ""I  _n2  “nk 


(6.3) 


Suppose  that  a  tracking  filter  is  applied  to  process  each  possible 
measurement  sequence.  The  a  posteriori  hypothesis  probability  of 
H  (k)  being  true  can  be  computed  recursively  using 

It 

p(z  (k)/H  (k-l),Zk"\ 


P (H  (k)/Z)  = 
mk 


p(z  (k)/Zk_1) 
"^k 


P  (H  (k-D/Z^1)  (6.4) 

mk-l 


k-1 

where  p(zn  (k)/Hj„  (k-D t  Z  )  is  the  probability  density  of  the 
k  k-1 
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residual  from  the  tracking  filter  using  H  (k-1)  and  z  (k)  . 

Hc-l  nk 

The  above  equation  can  be  derived  as  a  special  case  of  the  results 
presented  in  [22] ,  [23] .  The  final  set  of  tracks  (total  N)  can  be 
chosen  as  those  N  feasible  hypotheses  with  the  largest  hypothesis 
probabilities,  i.e., 

max  { P (H  (k)/Zk);  m  =  1,...,  Nk}  (6.5) 

mk  K 

{n;  h  (k  )c9} 
mk 

where  the  feasible  set#  9,  is  the  restriction  that  each  measure¬ 
ment  at  a  given  time  can  be  used  only  once,  i.e., 

9"  =  {H  (k)  :  H.  (k)  D  H.  (k)  =  <j>  for  i  ^  j}  (6.6) 

rak  k  ->k 

The  computational  requirement  of  the  above  method  is 
clearly  non-trivial.  In  fact,  the  above  optimization  problem 
defines  a  N-dimensional  assignment  problem.  To  the  best  of 
the  authors'  knowledge,  the  N-dimensional  extension  of  the  Hungari¬ 
an  algorithm  is  not  yet  available.  In  many  applications,  one 
may  be  able  to  pre-cluster  the  detections  so  that  the  search  over 
the  entire  set  of  detections  is  not  necessary.  Other  physical 
constraints  can  sometimes  be  imposed  to  reduce  the  search  require¬ 
ments  depending  upon  given  systems  and  application. 

A  similar  approach  using  a  maximum  likelihood  method 

*A  more  parametric  approach  for  modelling  this  probability  density 
function  is  given  in  Refs.  [40],  [42]  and  [46]  in  which  situations 
including  a  priori  traget  distribution  and  the  probability  of  a 
given  number  of  detections  were  also  considered. 
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was  described  in  [42]  in  which  the  multidimension  search 
problem  was  reduced  to  a  0-1  integer  programming  problem. 
Discussions  on  other  track  initiation  techniques  can  be  found 
in  [41],  [46  j  ,  [48]  . 

6.2  Track  Maintenance 

Once  track  files  have  been  established,  the  computation¬ 
al  requirement  is  greatly  reduced.  This  is  because  for  each 
track  file  one  is  only  required  to  search  the  "admissible" 
region  dictated  by  the  covariance  of  the  filter  residual  process. 

We  note  that  a  slightly  modified  method  of  the  track 
initiation  algorithm  discussed  in  6.1  can  be  applied  to  the 
track  maintenance  problem.  That  is,  one  establishes  a  new 
hypothesis  for  each  detection  resident  in  the  admissible  region. 
This  procedure  results  in  an  exponentially  growing  number  of 
track  files.  One  can  inhibit  the  growing  memory  and  computation¬ 
al  requirement  by  selecting  a  tree  depth  and  conducting  a 
global  search  for  a  set  of  feasible  tracks  having  the  highest 
hypothesis  probabilities  (eqs.  (6.5),  (6.6)).  Another  approach 
is  to  combine  a  set  of  "most  likely  hypothesis"  growing 
out  of  the  same  track  file  using  the  weighted  sum  of  state  esti¬ 
mates  with  the  hypothesis  probabilities  as  weighting  factors. 

This  second  approach  is  the  basis  of  the  Bayesian  tracker 
presented  by  Singer  et  al.[43j,  [44],  If  the  depth  is  equal  to 
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one,  i.e.,  one  combines  all  admissible  detections  at  each  scan, 
then  one  obtains  the  probabilistic  data  association  filter  of 
Bar-Sahlom  and  Tse  [45] .  We  emphasize  however,  that  the  approaches 
of  [43]- [45]  are  suitable  for  tracking  in  a  cluttered  environment 
(see  also  [47])  and  do  not  directly  address  the  multiple  target 
tracking  issue. 

The  above  discussion  did  not  include  problems  of  track 
termination,  imperfect  sensor  detection,  false  measurements,  etc. 
Reid  [46]  provides  an  interesting  paper  which  discusses  extensively 
features  of  both  track  initiation  and  track  maintenance. 

We  note  that  the  subject  of  tracking  in  a  dense  target 
environment  is  of  current  interest.  Although  all  quoted  refer¬ 
ences  present  general  approaches  to  this  problem  each  application 
typically  has  its  own  unique  features  which  impose  certain  res¬ 
trictions  resulting  in  a  variety  of  modifications.  Reference 
[48]  represents  such  a  situation  in  which  practical  constraints 
were  considered  for  tracking  algorithm  design  with  a  passive  sen- 


7.  ANALYSIS  TECHNIQUES 

In  system  design  and  trade-off  studies  involving 
target  tracking,  one  is  often  required  to  evaluate  the  performance 
of  trajectory  estimation  (tracking  and  prediction,  etc.)  for 
specific  situations.  Although  Monte  Carlo  simulation  is  a 
frequently  used  method,  it  is  nevertheless  time  consuming  and 
costly.  For  a  quick  and  general  performance  evaluation  par¬ 
ticularly  in  system  trade-off  studies  covariance  analysis  tech¬ 
niques  are  most  appropriate.  The  following  discussion  treats 
the  application  of  polynomial  analysis,  the  Riccati  Equation  and 
the  Cramer- Rao  bound. 

7.1  Polynomial  Analysis 

Let  Pq ,  vQ,  and  aQ  denote  at  t=0,  the  position,  velo¬ 
city  and  acceleration,  respectively,  of  a  moving  object  projected 
along  a  given  coordinate  (it  may  be  either  range  or  angles) 
then  its  position  at  an  arbitary  time  t  is 

p(t)  =  pQ  +  Vq t  +  j  aQt2.  (7.1) 

Let  p(tn) ,  n=l,...,  N  denote  a  set  of  noisy  measurements  of 
p(tn).  For  the  purpose  of  convenience,  let  the  time  samples  be 
taken  uniformly  with  spacing  T  and  let  the  total  time  interval  be 
centered  about  t=0.  Under  these  conditions,  the  estimates  of 
target  position,  velocity,  and  acceleration  at  an  arbitrary  time 
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t  are  given  by 
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The  covariance  of  the  state  estimate  (Eq.  (7.2))  is 


P  =  a2<S>( t)S  $(t)T  (7.6) 

2 

where  a  is  the  variance  of  the  measurement  noise. 

We  note  that  the  above  result,  especially  (7.6),  is 
very  useful  because  of  its  simplicity  and  applicability  to  a  wide 
range  of  problems.  Specific  application  of  the  above  analysis  to 
re-entry  vehicle  tracking  problems  and  discussions  on  polynomial 
analyses  are  found  in  [49] . 

7.2  The  Use  of  the  Riccati  Equation 

The  above  analysis  is  straightforward  and  its  cal¬ 
culation  can  be  carried  out  using  desk  top  (or  pocket) 
calculators.  The  drawback  is  that  it  becomes  overly  optimistic 
for  long  track  intervals  ((N-l)T).  This  results  since  the 
target  dynamic  model  (Eq.  (7.1))  does  not  include  process  noise. 
To  circumvent  this  problem,  one  may  compute  the  covariance 
matrix  using  the  Riccati  equation  (Eq.  (2.13)).  This 
will  require  a  computer  when  nonlinear  dynamics  and/or  measure¬ 
ment  equations  are  involved.  It  is  nevertheless  a  convenient 
method  since  the  Monte  Carlo  simulation  is  not  required. 

We  note  that  the  use  of  Riccati  equation  for  nonlinear 
problems  only  represents  an  approximate  error  analysis.  If  one 
choses  a  sufficiently  large  process  noise  covariance,  then  the 
result  of  the  Riccati  equation  is  an  upper  bound  to  the  actual 
performance,  [9],  The  tightness  of  this  bound  is  however,  rather 
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difficult  to  assess 


7. 3  The  Use  of  The  Cramer-Rao  Lower  Bound 


There  exist  some  situations  in  which  the  process  noise 
is  negligible.  The  techniques  of  Section  7.1  are  therefore 
applicable  to  a  limited  degree  although  they  do  ignore  the  coupl¬ 
ing  between  coordinates  and  assume  a  linear  measurement  system. 

When  the  coupling  between  coordinates  becomes  significant  one 
must  use  the  Riccati  equation  to  compute  the  covariance.  When 
the  process  noise  term  is  indeed  negligible  (e.g. ,  the  exo- 
atmospheric  trajectory  estimation,  Section  5.1),  the  solution  of 
the  Riccati  equation  using  Jacobian  matrices  evaluated  along  the 
true  trajectory  becomes  the  Cramer-Rao  lower  bound  on  the  covariance 
of  the  state  estimates,  [49],  [50],  [51],  Furthermore,  the  Riccati 
equation  has  a  closed  form  solution.  As  an  example,  for  the  dis¬ 
crete  case. 


r  n 

?  =  I  Zrf  F  TH  TR"1H  F 

N  I  *■'*«  n  n  n  n  n 
L  n=l 


-1 


(7.7) 


N-l 


rn  *  n  V1 


for  n  <  N-l 


n=l 


(7.8) 


=  I 


for  n  =  N 


where  PN  is  the  Cramer-Rao  bound  for  the  covariance  of  Xj^,  1  is 
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am  identity  matrix  and  $n  and  Hn  are  system  (discrete)  and  measure¬ 
ment  Jacobian  matrices  evaluated  at  xn,  respectively. 

The  significance  of  the  above  fact  is  that  the  actual 
performance  is  lower  bounded  by  the  solution  of  the  Riccati 
equation.  An  extended  Kalman  filter  generally  does  not  perform 
as  well  while  a  properly  designed  filter  can  asymptotically 
achieve  this  bound,  see  [30],  [57]. 

Finally  we  remark  that  the  corresponding  Cramer-Rao 
bound  for  systems  with  process  noise  is  very  difficult  to  compute. 
Versions  of  lower  bounds  to  the  Cramer-Rao  bound  have  been  pro¬ 
posed  [51],  [52],  these  bounds  are  tight  only  for  very  large  signal 
to  noise  ratios. 

8.  CONCLUDING  REMARKS 


In  this  report,  we  have  presented  a  survey  of  problems 
and  solutions  which  deal  with  target  tracking.  Although  this 
problem  has  been  of  active  concern  to  practitioners  in  both 
military  and  civilian  applications  for  many  years,  new  problems 
still  emerge  and  challenge  systems  engineers.  It  is  felt  that  the 
basic  approach  to  tracking  filter  design  is  no  longer  an  issue. 
However,  the  problem  of  integrating  the  tracking  system  into  the 
overall  command,  control,  and  communication  structure  to  achieve 
improved  performance  while  minimizing  data  processing  requirements 
represent  the  nature  of  problems  that  require  current  attention. 
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